#the following script allows one to replicate all weekly analyses pertaining to Cabo Delgado
#it uses the following files: weekly_cyclone_ts.dta and pests.r

#load necessary libraries
library(foreign)
library(MASS)
library(astsa)

#source in pests R code
source("pests.r")

#set random number seed
set.seed(10)

#read-in dataset
cyclonedata <- read.dta(file = "weekly_cyclone_ts.dta")

#attach data
attach(cyclonedata)

#summary statistics for table A.1
summary(moz_rebcit_acled)
sd(moz_rebcit_acled,na.rm=TRUE)
summary(fullcyclone)
sd(fullcyclone,na.rm=TRUE)
summary(moz_govreb_acled)
sd(moz_govreb_acled,na.rm=TRUE)
summary(lagmoz_govreb_acled)
sd(lagmoz_govreb_acled,na.rm=TRUE)
summary(moz_rebgov_acled)
sd(moz_rebgov_acled,na.rm=TRUE)
summary(lagmoz_rebgov_acled)
sd(lagmoz_rebgov_acled,na.rm=TRUE)
summary(lag2moz_rebgov_acled)
sd(lag2moz_rebgov_acled,na.rm=TRUE)
summary(moz_govcit_acled)
sd(moz_govcit_acled,na.rm=TRUE)
summary(lagmoz_govcit_acled)
sd(lagmoz_govcit_acled,na.rm=TRUE)
summary(moz_ncd)
sd(moz_ncd,na.rm=TRUE)
summary(lagmoz_ncd)
sd(lagmoz_ncd,na.rm=TRUE)
summary(lag2moz_ncd)
sd(lag2moz_ncd,na.rm=TRUE)
summary(russian)
sd(russian,na.rm=TRUE)

##################
##Appendix Plots##
##################

#Figure A.3a-A.3b
par(mfrow=c(1,1))
plot(moz_rebcit_acled, ylab="Count", main="Weekly Counts")
acftmp <- counts.acf(moz_rebcit_acled)
tsplot(moz_rebcit_acled)
tsplot(moz_rebcit_acled, xaxt="n",ylab="Insurgent Violence Against Civilians",xlab="Weekly Observations: March 2018 - April 2020")
axis(1, at=c(0, 20, 40, 60, 80, 100),labels=c("02/2018", "07/2018", "12/2018", "04/2019", "09/2019", "01/2020")) 

#Figure A.4a-A.4b
par(mfrow=c(1,1))
acftmp <- counts.acf(moz_rebcit_acled[1:54])
tsplot(moz_rebcit_acled[1:54], xaxt="n",ylab="Insurgent Violence Against Civilians",xlab="Weekly Observations: March 2018 - March 2019")
axis(1, at=c(0, 10, 20, 30, 40, 50),labels=c("02/2018", "05/2018", "07/2018", "09/2018", "10/2018", "02/2019")) 

#Figure A.5a-A.4b
par(mfrow=c(1,1))
plot(moz_rebcit_acled[55:114], ylab="Count", main="Weekly Counts")
acftmp <- counts.acf(moz_rebcit_acled[55:114])
tsplot(moz_rebcit_acled[55:114], xaxt="n",ylab="Insurgent Violence Against Civilians",xlab="Weekly Observations: March 2019 - April 2020")
axis(1, at=c(0, 10, 20, 30, 40, 50,60),labels=c("03/2019", "05/2019", "08/2019", "10/2019", "12/2019", "03/2020", "04/2020")) 


############
##Select p##
############

#for moz_rebcit_acled
# PAR(1) Model
PAR1 <- Parp(moz_rebcit_acled ~ 1, p=1)
print(PAR1)
# PAR(2) Model
PAR2 <- Parp(moz_rebcit_acled ~ 1, p=2)
print(PAR2)
# PAR(3) Model
PAR3 <- Parp(moz_rebcit_acled ~ 1, p=3)
print(PAR3)
# PAR(4) Model
PAR4 <- Parp(moz_rebcit_acled ~ 1, p=4)
print(PAR4)
# PAR(5) Model
PAR5 <- Parp(moz_rebcit_acled ~ 1, p=5)
print(PAR5)
# PAR(6) Model
PAR6<- Parp(moz_rebcit_acled ~ 1, p=6)
print(PAR6)


#########################################
##Select proper exogenous variable lags##
#########################################

#small
Cyclone.PAR2 <- Parp(moz_rebcit_acled ~ fullcyclone+moz_govreb_acled, p=2)
print(Cyclone.PAR2)
Cyclone.PAR2 <- Parp(moz_rebcit_acled ~ fullcyclone+moz_govreb_acled+lagmoz_govreb_acled, p=2)
print(Cyclone.PAR2)
Cyclone.PAR2 <- Parp(moz_rebcit_acled ~ fullcyclone+moz_govreb_acled+lagmoz_govreb_acled+lag2moz_govreb_acled, p=2)
print(Cyclone.PAR2)
Cyclone.PAR2 <- Parp(moz_rebcit_acled ~ fullcyclone+lagfullcyclone+moz_govreb_acled, p=2)
print(Cyclone.PAR2)
Cyclone.PAR2 <- Parp(moz_rebcit_acled ~ fullcyclone+lagfullcyclone+moz_govreb_acled+lagmoz_govreb_acled, p=2)
print(Cyclone.PAR2)
Cyclone.PAR2 <- Parp(moz_rebcit_acled ~ fullcyclone+lagfullcyclone+moz_govreb_acled+lagmoz_govreb_acled+lag2moz_govreb_acled, p=2)
print(Cyclone.PAR2)
Cyclone.PAR2 <- Parp(moz_rebcit_acled ~ fullcyclone+lagfullcyclone+lag2fullcyclone+moz_govreb_acled, p=2)
print(Cyclone.PAR2)
Cyclone.PAR2 <- Parp(moz_rebcit_acled ~ fullcyclone+lagfullcyclone+lag2fullcyclone+moz_govreb_acled+lagmoz_govreb_acled, p=2)
print(Cyclone.PAR2)
Cyclone.PAR2 <- Parp(moz_rebcit_acled ~ fullcyclone+lagfullcyclone+lag2fullcyclone+moz_govreb_acled+lagmoz_govreb_acled+lag2moz_govreb_acled, p=2)
print(Cyclone.PAR2)

#medium
Cyclone.PAR2 <- Parp(moz_rebcit_acled ~ fullcyclone+moz_govreb_acled+lagmoz_govreb_acled+moz_rebgov_acled, p=2)
print(Cyclone.PAR2)
Cyclone.PAR2 <- Parp(moz_rebcit_acled ~ fullcyclone+moz_govreb_acled+lagmoz_govreb_acled+moz_rebgov_acled+lagmoz_rebgov_acled, p=2)
print(Cyclone.PAR2)
Cyclone.PAR2 <- Parp(moz_rebcit_acled ~ fullcyclone+moz_govreb_acled+lagmoz_govreb_acled+moz_rebgov_acled+lagmoz_rebgov_acled+lag2moz_rebgov_acled, p=2)
print(Cyclone.PAR2)

#large
Cyclone.PAR2 <- Parp(moz_rebcit_acled ~ fullcyclone+moz_govreb_acled+lagmoz_govreb_acled+moz_rebgov_acled+lagmoz_rebgov_acled+lag2moz_rebgov_acled+moz_govcit_acled, p=2)
print(Cyclone.PAR2)
Cyclone.PAR2 <- Parp(moz_rebcit_acled ~ fullcyclone+moz_govreb_acled+lagmoz_govreb_acled+moz_rebgov_acled+lagmoz_rebgov_acled+lag2moz_rebgov_acled+moz_govcit_acled+lagmoz_govcit_acled, p=2)
print(Cyclone.PAR2)
Cyclone.PAR2 <- Parp(moz_rebcit_acled ~ fullcyclone+moz_govreb_acled+lagmoz_govreb_acled+moz_rebgov_acled+lagmoz_rebgov_acled+lag2moz_rebgov_acled+moz_govcit_acled+lagmoz_govcit_acled+lag2moz_govcit_acled, p=2)
print(Cyclone.PAR2)

#extra large
Cyclone.PAR2 <- Parp(moz_rebcit_acled ~ fullcyclone+moz_govreb_acled+lagmoz_govreb_acled+moz_rebgov_acled+lagmoz_rebgov_acled+lag2moz_rebgov_acled+moz_govcit_acled+lagmoz_govcit_acled+moz_ncd+russian, p=2)
print(Cyclone.PAR2)
Cyclone.PAR2 <- Parp(moz_rebcit_acled ~ fullcyclone+moz_govreb_acled+lagmoz_govreb_acled+moz_rebgov_acled+lagmoz_rebgov_acled+lag2moz_rebgov_acled+moz_govcit_acled+lagmoz_govcit_acled+moz_ncd+lagmoz_ncd+russian, p=2)
print(Cyclone.PAR2)
Cyclone.PAR2 <- Parp(moz_rebcit_acled ~ fullcyclone+moz_govreb_acled+lagmoz_govreb_acled+moz_rebgov_acled+lagmoz_rebgov_acled+lag2moz_rebgov_acled+moz_govcit_acled+lagmoz_govcit_acled+moz_ncd+lagmoz_ncd+lag2moz_ncd+russian, p=2)
print(Cyclone.PAR2)
Cyclone.PAR2 <- Parp(moz_rebcit_acled ~ fullcyclone+moz_govreb_acled+lagmoz_govreb_acled+moz_rebgov_acled+lagmoz_rebgov_acled+lag2moz_rebgov_acled+moz_govcit_acled+lagmoz_govcit_acled+moz_ncd+russian+lagrussian, p=2)
print(Cyclone.PAR2)
Cyclone.PAR2 <- Parp(moz_rebcit_acled ~ fullcyclone+moz_govreb_acled+lagmoz_govreb_acled+moz_rebgov_acled+lagmoz_rebgov_acled+lag2moz_rebgov_acled+moz_govcit_acled+lagmoz_govcit_acled+moz_ncd+lagmoz_ncd+russian+lagrussian, p=2)
print(Cyclone.PAR2)
Cyclone.PAR2 <- Parp(moz_rebcit_acled ~ fullcyclone+moz_govreb_acled+lagmoz_govreb_acled+moz_rebgov_acled+lagmoz_rebgov_acled+lag2moz_rebgov_acled+moz_govcit_acled+lagmoz_govcit_acled+moz_ncd+lagmoz_ncd+lag2moz_ncd+russian+lagrussian, p=2)
print(Cyclone.PAR2)
Cyclone.PAR2 <- Parp(moz_rebcit_acled ~ fullcyclone+moz_govreb_acled+lagmoz_govreb_acled+moz_rebgov_acled+lagmoz_rebgov_acled+lag2moz_rebgov_acled+moz_govcit_acled+lagmoz_govcit_acled+moz_ncd+russian+lagrussian+lag2russian, p=2)
print(Cyclone.PAR2)
Cyclone.PAR2 <- Parp(moz_rebcit_acled ~ fullcyclone+moz_govreb_acled+lagmoz_govreb_acled+moz_rebgov_acled+lagmoz_rebgov_acled+lag2moz_rebgov_acled+moz_govcit_acled+lagmoz_govcit_acled+moz_ncd+lagmoz_ncd+russian+lagrussian+lag2russian, p=2)
print(Cyclone.PAR2)
Cyclone.PAR2 <- Parp(moz_rebcit_acled ~ fullcyclone+moz_govreb_acled+lagmoz_govreb_acled+moz_rebgov_acled+lagmoz_rebgov_acled+lag2moz_rebgov_acled+moz_govcit_acled+lagmoz_govcit_acled+moz_ncd+lagmoz_ncd+lag2moz_ncd+russian+lagrussian+lag2russian, p=2)
print(Cyclone.PAR2)


###########
##Table 1##
###########

#small specification
Cyclone.small <- Parp(moz_rebcit_acled ~ fullcyclone+moz_govreb_acled+lagmoz_govreb_acled, p=2)
print(Cyclone.small)

#medium specification
Cyclone.medium <- Parp(moz_rebcit_acled ~ fullcyclone+moz_govreb_acled+lagmoz_govreb_acled+moz_rebgov_acled+lagmoz_rebgov_acled+lag2moz_rebgov_acled, p=2)
print(Cyclone.medium)

#large specification
Cyclone.large <- Parp(moz_rebcit_acled ~ fullcyclone+moz_govreb_acled+lagmoz_govreb_acled+moz_rebgov_acled+lagmoz_rebgov_acled+lag2moz_rebgov_acled+moz_govcit_acled+lagmoz_govcit_acled, p=2)
print(Cyclone.large)

#extra large specification
Cyclone.xl <- Parp(moz_rebcit_acled ~ fullcyclone+moz_govreb_acled+lagmoz_govreb_acled+moz_rebgov_acled+lagmoz_rebgov_acled+lag2moz_rebgov_acled+moz_govcit_acled+lagmoz_govcit_acled+moz_ncd+lagmoz_ncd+lag2moz_ncd+russian, p=2)
print(Cyclone.xl)

###########
##Table 2##
###########

# Estimation of the dynamic multipliers for one unit changes in each
# covariate.  These compute the effects of changes in each variable
# holding the others at their means.

small.mults <- montecarlo.parp.multipliers(Cyclone.small,
                                             xvec=rep(1,4), n=1000)

medium.mults <- montecarlo.parp.multipliers(Cyclone.medium,
                                             xvec=rep(1,7), n=1000)

large.mults <- montecarlo.parp.multipliers(Cyclone.large, xvec=rep(1,9),
                                             n=1000)

xl.mults <- montecarlo.parp.multipliers(Cyclone.xl, xvec=rep(1,13),
                                             n=1000)

# To summarize the simulated multipliers, you need to sum over the
# simulated values.  Here we summarize the 0.16, 0.5 and 0.84
# quantiles.  The columns of this output correspond to each of the
# covariates (with the intercept as the first).  The first matrix are
# the short run or impact multipliers, the second are the long run or
# total multipliers

apply(small.mults, c(2,3), quantile, probs=c(0.16, 0.5, 0.84))
apply(medium.mults, c(2,3), quantile, probs=c(0.16, 0.5, 0.84))
apply(large.mults, c(2,3), quantile, probs=c(0.16, 0.5, 0.84))
apply(xl.mults, c(2,3), quantile, probs=c(0.16, 0.5, 0.84))


#############
##Table A.2##
#############

#small specification
Cyclone.PAR2 <- Parp(moz_rebcit_acled ~ cyclonekenneth+moz_govreb_acled+lagmoz_govreb_acled, p=2)
print(Cyclone.PAR2)

#medium specification
Cyclone.PAR2 <- Parp(moz_rebcit_acled ~ cyclonekenneth+moz_govreb_acled+lagmoz_govreb_acled+moz_rebgov_acled+lagmoz_rebgov_acled+lag2moz_rebgov_acled, p=2)
print(Cyclone.PAR2)

#large specification
Cyclone.PAR2 <- Parp(moz_rebcit_acled ~ cyclonekenneth+moz_govreb_acled+lagmoz_govreb_acled+moz_rebgov_acled+lagmoz_rebgov_acled+lag2moz_rebgov_acled+moz_govcit_acled+lagmoz_govcit_acled, p=2)
print(Cyclone.PAR2)

#extra large specification
Cyclone.PAR2<- Parp(moz_rebcit_acled ~ cyclonekenneth+moz_govreb_acled+lagmoz_govreb_acled+moz_rebgov_acled+lagmoz_rebgov_acled+lag2moz_rebgov_acled+moz_govcit_acled+lagmoz_govcit_acled+moz_ncd+lagmoz_ncd+lag2moz_ncd+russian, p=2)
print(Cyclone.PAR2)


#############
##Table A.3##
#############

Cyclone.PAR2 <- Parp(moz_rebcit_acled ~ fullcyclone+shortcyclone+moz_govreb_acled+lagmoz_govreb_acled, p=2)
print(Cyclone.PAR2)

#medium specification
Cyclone.PAR2 <- Parp(moz_rebcit_acled ~ fullcyclone+shortcyclone+moz_govreb_acled+lagmoz_govreb_acled+moz_rebgov_acled+lagmoz_rebgov_acled+lag2moz_rebgov_acled, p=2)
print(Cyclone.PAR2)

#large specification
Cyclone.PAR2 <- Parp(moz_rebcit_acled ~ fullcyclone+shortcyclone+moz_govreb_acled+lagmoz_govreb_acled+moz_rebgov_acled+lagmoz_rebgov_acled+lag2moz_rebgov_acled+moz_govcit_acled+lagmoz_govcit_acled, p=2)
print(Cyclone.PAR2)

#extra large specification
Cyclone.PAR2<- Parp(moz_rebcit_acled ~ fullcyclone+shortcyclone+moz_govreb_acled+lagmoz_govreb_acled+moz_rebgov_acled+lagmoz_rebgov_acled+lag2moz_rebgov_acled+moz_govcit_acled+lagmoz_govcit_acled+moz_ncd+lagmoz_ncd+lag2moz_ncd+russian, p=2)
print(Cyclone.PAR2)

#############
##Table A.4##
#############

#small specification
Cyclone.PAR2 <- Parp(moz_rebcit_acled ~ fullcyclone+moz_govreb_acled+lagmoz_govreb_acled, p=1)
print(Cyclone.PAR2)

#medium specification
Cyclone.PAR2 <- Parp(moz_rebcit_acled ~ fullcyclone+moz_govreb_acled+lagmoz_govreb_acled+moz_rebgov_acled+lagmoz_rebgov_acled+lag2moz_rebgov_acled, p=1)
print(Cyclone.PAR2)

#large specification
Cyclone.PAR2 <- Parp(moz_rebcit_acled ~ fullcyclone+moz_govreb_acled+lagmoz_govreb_acled+moz_rebgov_acled+lagmoz_rebgov_acled+lag2moz_rebgov_acled+moz_govcit_acled+lagmoz_govcit_acled, p=1)
print(Cyclone.PAR2)

Cyclone.PAR2 <- Parp(moz_rebcit_acled ~ fullcyclone+moz_govreb_acled+lagmoz_govreb_acled+moz_rebgov_acled+lagmoz_rebgov_acled+lag2moz_rebgov_acled+moz_govcit_acled+lagmoz_govcit_acled+moz_ncd+lagmoz_ncd+lag2moz_ncd+russian, p=1)
print(Cyclone.PAR2)

#############
##Table A.5##
#############

#small specification
Cyclone.PAR2 <- Parp(moz_rebcit_acled ~ fullcyclone+moz_govreb_acled+lagmoz_govreb_acled, p=3)
print(Cyclone.PAR2)

#medium specification
Cyclone.PAR2 <- Parp(moz_rebcit_acled ~ fullcyclone+moz_govreb_acled+lagmoz_govreb_acled+moz_rebgov_acled+lagmoz_rebgov_acled+lag2moz_rebgov_acled, p=3)
print(Cyclone.PAR2)

#large specification
Cyclone.PAR2 <- Parp(moz_rebcit_acled ~ fullcyclone+moz_govreb_acled+lagmoz_govreb_acled+moz_rebgov_acled+lagmoz_rebgov_acled+lag2moz_rebgov_acled+moz_govcit_acled+lagmoz_govcit_acled, p=3)
print(Cyclone.PAR2)

Cyclone.PAR2 <- Parp(moz_rebcit_acled ~ fullcyclone+moz_govreb_acled+lagmoz_govreb_acled+moz_rebgov_acled+lagmoz_rebgov_acled+lag2moz_rebgov_acled+moz_govcit_acled+lagmoz_govcit_acled+moz_ncd+lagmoz_ncd+lag2moz_ncd+russian, p=3)
print(Cyclone.PAR2)

#############
##Table A.6##
#############

#small specification
Cyclone.PAR2 <- Parp(moz_rebcit_acled ~ fullcyclone+moz_govreb_acled+lagmoz_govreb_acled+is, p=2)
print(Cyclone.PAR2)

#medium specification
Cyclone.PAR2 <- Parp(moz_rebcit_acled ~ fullcyclone+moz_govreb_acled+lagmoz_govreb_acled+moz_rebgov_acled+lagmoz_rebgov_acled+lag2moz_rebgov_acled+is, p=2)
print(Cyclone.PAR2)

#large specification
Cyclone.PAR2 <- Parp(moz_rebcit_acled ~ fullcyclone+moz_govreb_acled+lagmoz_govreb_acled+moz_rebgov_acled+lagmoz_rebgov_acled+lag2moz_rebgov_acled+moz_govcit_acled+lagmoz_govcit_acled+is, p=2)
print(Cyclone.PAR2)

Cyclone.PAR2 <- Parp(moz_rebcit_acled ~ fullcyclone+moz_govreb_acled+lagmoz_govreb_acled+moz_rebgov_acled+lagmoz_rebgov_acled+lag2moz_rebgov_acled+moz_govcit_acled+lagmoz_govcit_acled+moz_ncd+lagmoz_ncd+lag2moz_ncd+russian+is, p=2)
print(Cyclone.PAR2)
 

#############
##Table A.7##
#############

#small specification
m1<- glm(moz_rebcit_acled ~ fullcyclone+moz_govreb_acled+lagmoz_govreb_acled, family="poisson")
summary(m1) 
logLik(m1)

#medium specification
m2<- glm(moz_rebcit_acled ~ fullcyclone+moz_govreb_acled+lagmoz_govreb_acled+moz_rebgov_acled+lagmoz_rebgov_acled+lag2moz_rebgov_acled, family="poisson")
summary(m2) 
logLik(m2)

#large specification
m3 <- glm(moz_rebcit_acled ~ fullcyclone+moz_govreb_acled+lagmoz_govreb_acled+moz_rebgov_acled+lagmoz_rebgov_acled+lag2moz_rebgov_acled+moz_govcit_acled+lagmoz_govcit_acled, family="poisson")
summary(m3)
logLik(m3)

#large specification
m4 <- glm(moz_rebcit_acled ~ fullcyclone+moz_govreb_acled+lagmoz_govreb_acled+moz_rebgov_acled+lagmoz_rebgov_acled+lag2moz_rebgov_acled+moz_govcit_acled+lagmoz_govcit_acled+moz_ncd+lagmoz_ncd+lag2moz_ncd+russian, family="poisson")
summary(m4)
logLik(m4)

#############
##Table A.8##
#############

#small specification
m1<- glm.nb(moz_rebcit_acled ~ fullcyclone+moz_govreb_acled+lagmoz_govreb_acled)
summary(m1) 
logLik(m1)

#medium specification
m2<- glm.nb(moz_rebcit_acled ~ fullcyclone+moz_govreb_acled+lagmoz_govreb_acled+moz_rebgov_acled+lagmoz_rebgov_acled+lag2moz_rebgov_acled)
summary(m2) 
logLik(m2)

#large specification
m3<- glm.nb(moz_rebcit_acled ~ fullcyclone+moz_govreb_acled+lagmoz_govreb_acled+moz_rebgov_acled+lagmoz_rebgov_acled+lag2moz_rebgov_acled+moz_govcit_acled+lagmoz_govcit_acled)
summary(m3)
logLik(m3)

#large specification
m4<- glm.nb(moz_rebcit_acled ~ fullcyclone+moz_govreb_acled+lagmoz_govreb_acled+moz_rebgov_acled+lagmoz_rebgov_acled+lag2moz_rebgov_acled+moz_govcit_acled+lagmoz_govcit_acled+moz_ncd+lagmoz_ncd+lag2moz_ncd+russian)
summary(m4)
logLik(m4)

##############
##Table A.10##
##############

#small specification (main setup)
Cyclone.small <- Parp(moz_rebcit_acled ~ fullcyclone, p=2)
print(Cyclone.small)

#small specification, kenneth variable
Cyclone.PAR2 <- Parp(moz_rebcit_acled ~ cyclonekenneth, p=2)
print(Cyclone.PAR2)

#small specification, short cyclone
Cyclone.PAR2 <- Parp(moz_rebcit_acled ~ fullcyclone+shortcyclone, p=2)
print(Cyclone.PAR2)

#small specification (alt p)
Cyclone.PAR2 <- Parp(moz_rebcit_acled ~ fullcyclone, p=1)
print(Cyclone.PAR2)

#small specification (alt p 2)
Cyclone.PAR2 <- Parp(moz_rebcit_acled ~ fullcyclone, p=3)
print(Cyclone.PAR2)

#small specification
Cyclone.PAR2 <- Parp(moz_rebcit_acled ~ fullcyclone+is, p=2)
print(Cyclone.PAR2)

#small specification, poisson
m1<- glm(moz_rebcit_acled ~ fullcyclone, family="poisson")
summary(m1) 
logLik(m1)

#small specification, negbin
m1<- glm.nb(moz_rebcit_acled ~ fullcyclone)
summary(m1) 
logLik(m1)

